UNCLASSIFIED 


AD  NUMBER 


AD886376 


NEW  LIMITATION  CHANGE 
TO 

Approved  for  public  release,  distribution 
unlimited 


FROM 

Distribution  authorized  to  U.S.  Gov't, 
agencies  and  their  contractors;  Critical 
Technology;  30  JUN  1971.  Other  requests 
shall  be  referred  to  Air  Force  Technical 
Applications  Center,  Patrick  AFB,  FL. 


AUTHORITY 


AFTAC  ltr ,  1  May  1972 


THIS  PAGE  IS  UNCLASSIFIED 


AD886376 


Computation  of  Multipie-Event  Probabilities, 


TELEDYNE  GEOTECH  ALEXANDRIA  VA 


30  JUN  1971 


COMPUTATION  OF  MULTIPLE-EVENT  PROBABILITIES 


SEISMIC  DATA  LABORATORY  REPORT  NO.  277 


AFTAC  Proj  ect  No . : 
Project  Title: 

ARPA  Order  No.: 

ARPA  Program  Code  No. 


Name  of  Contractor: 


Contract  No. : 

Date  of  Contract: 

Amount  of  Contract: 
Contract  Expiration  Date: 
Project  Manager: 


VELA  T/ 2706 

Seismic  Data  Laboratory 


2F-10 


TELEDYNE  GEOTECH 


F33657-72-C-0009 

01  July  1971 

$  1,290,000 

30  June  1972 

Royal  A.  Hartenberger 
(7031  836-7647 


P.  u.  Box  334,  Alexandria,  Virginia 


This  document  is  subject  to  special  export  controls  and  each 
transmittal  to  foreign  governments  or  foreign  nationals  may 
be  made  only  with  prior  approval  of  Chief,  AFTAC. 


f  t  ^  t-.  >**'>' 


VWffK  « 


Unc lass  if  if  d 


Security  Classification 


DOCUMENT  CONTROL  DATA  •  R&O 

(Security  c/««»//fCJ</on  of  i/ffo.  C*>dy  of  obofrocf  ond  fnd*«in|  annofaffon  muil  bo  anfarad  w#>«n  W»a  ovaraJJ  rtporf  i«  cJaaatffadj 


I  OpTgIMATINC  ACTIVITY  fCorpor.l.  author] 

-  TELEDYNE  GE0TEC11 
.ALEXANDRIA,  VIRGINIA 


2a  pkpopt  »icupitv  c  las»ipic» tion 

Unclass i f ied 


26  anoup 


1  PEPOPT  TITLE 


COMPUTATION  OF  MULTIPLE- EVENT  PROBABILITIES 


a  oeicriptive  n6tcs  TTrpa  a  niSS  mi  mcSmTw  <(•<••) 

Scientific 


:  AUT^CS,'!,’ fir*' naaia.  fl-l  inftltl) 

Wirth,  Mark  H. 


S.  PEPOPT  DATE 


6/3U/71 


7«>  TOTAL  NO.  OP  P Adit 


14 


76.  Np  OP  PSP1 

3 


Sa  CONTRACT  OP  4PANT  HP. 


nFJUsxrj*c**f 

k.  PPOJSCT  NO. 

VELA  T/0706 
eARPA  Order  No.  024 

<*ARPA  Program  Code  No.  9F10 


Sa  OPiaiNATOP'S  PtPOPT  NUMaaPCJJ 


277 


tb  OTHER  REPORT  MOffJ  (A  vy  otbar  numb«r«  that  may  ba  aaafdnad 
di /•  rmpart) 


10  A  V  A  1L  A1ILITY/LIMITATIOM  NOTICE* 

This  aocument  is  subject  to  special  export  controls  and  each  trans 
mittal  to  foreign  governments  or  foreign  nationals  may  be  made  onl)| 
with  prior  approval  of  Chief-  AFT AC 


11  SUPPLEMENTAPV  NOTES 


12  SPONSOPINO  MILITARY  ACTIVITY 

Advanced  Research  Projects  Agency 
Nuclear  Monitoring  Research  Office 

■Waslnaglan...  P.t— _ 


tl  ABSTRACT 


A  new,  simple  and  efficient  algorithm  is  described 
which  computes  the  probability  of  occurrence  of  at  least 
(or  exactly)  k  out  of  N  events,  given  the  individual  event 
probabilities.  A  comparison  is  made  with  two  previous 
methous.  Computer  program  listings  are  given. 
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ABSTRACT 

A  new,  simple  and  efficient  algorithm  is  described 
which  computes  the  probability  of  occurrence  of  at 
least  (or  exactly)  k/ out  of  N  events,  given  the  indi¬ 
vidual  event  probabilities.  A  comparison  is  made  with 
two  previous  methods.  Computer  program  listings  are 
given. 
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TilC  ALOUKliilM 


It  is  desired  to  compute  the  probability  of  at 
least  k  out  of  N  events  occurring,  given  the  individual 
event  probabilities  p^.  This  can  be  obtained  by  summa¬ 
tion  from  the  probability  of  exactly  k  events,  P(k), 

i.e. 

N 

F(=k)  =  2]  PU) 
i=k 


Thus  we  need  only  compute  the  probability  of  exactly  k 
events . 

The  algorithm  is  based  on  the  observation  that  given 
two  groups  of  events,  with  their  associated  probabilities 
of  occurrence,  the  probability  of  exactly  k  total  out  of 
the  two.  groups  can  be  found  by 


N 

where  by  P(jJS)  we  denote  the  probability  of  exactly  k 
out  of  N  events  from  the  set  S.  The  index  i  is  subject 
to  the  limitations  0  <_  i  <  and  0  £  k-i  <_  N2,  i.e. 
sum  over  i  -=  max(0  ,k-N2)  , . . .  .minCk,?^) .  Note  that 
PCj|S  )  =  Pj  ,  and  •  P(j!j  |  Sj )  =  1-p^,  by  definition.  Thus 
we  are  initially  given  the  probabilities  for  groups  of 
size  1.  Using  the  combining  rule  (2),  we  can  recursively 
compute  the  probabilities  for  groups  of  size  2,4,8, ...,N. 
Since-*C2)  work?  for  disparate  group  sizes,  no  problem 
is  caused  by  the  odd  groups  that  result  when  N  is  not  a 


power  of  t“o;  An  odd  group  without  a  mate  to  pa? r  with 
at  any  stage  is  simply  passed  on  unchanged  to  the  next 
stage:  at  a  later  stage  it  will  be  paired  with  a  group 
haying  a  different  number  of  elements.  A  schematic 
illustration  of  the  operation  of  the  .algorithm- for 
N  =  5:  1 


input: 


output: 


11111 

0101010101 

012-012-01 

01234—01 

012345 - 

12345 


where  the  numerals  indicate  the  index  of  and  the 

positions  indicate  the  relative  positions  in  an  array 
of  storage,  i.e.  each  "2"  in  the  diagram  represents  the 
location  of  a  probability  of  exactly  two  events,  etc. 
Dashes  indicate  blanks  (unused  locations).  The  final 
summation  is  not  shown. 

A  listing  of  a  computer  subroutine  (NETPROB  III) 
employing  the  algorithm  is  given  in  the  Appendix.  The 
restriction  N  <_  50  is  completely  arbitrary  for  this 
routine  and  derives  from  the  requirement  of  at  least 
4N  locations  in  the  working  array.  The  routine  can  be 
modified  to  output  the  probability  of  exactly  k  events 
by  eliminating  the  final  accumulation  (statement  200). 
Note  that  the  P(>_k)  are  computed  for  all  k  at  once  by 
this  algorithm  (k  =  1,...,N). 
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TIMING 


It  is  easily  shown  that  the  combining  ooeration  (2) 
2 

requires  ~M  /4  operations  for  M  =  N^  +  =  2,4,8,...,N. 

Thus  the  total  time  should  be  proportional  to 


4.N  16  N  64. N 

4  2  +  4*4  +  4  8  + 


=  |  ( 1+2+4+ •  *  *+~) 


mzii  „  n2 


Timing  tests  on  a  CDC  1604B  c  mputer  gave  a  timing  curve 


T  «  .08  (N+25)N  msec. 


h> 
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REMARKS 

2  * 

Any  N  process  for  which  a  combining  rule  analogous 

s 

to  (2)  can  be  given  requiring  on^y  M  operations,  where  M 

is  the  total  number^of  elements  in  the  two  groups,  can  bc^ 

converted  to  an  N  log2N  process  by  this  -binary  subdivi-'- 

sion  algorithm..  In  fact,  the  author  has  used  a  virtually 

identical  subdivision  routine  Sas  the  outer  routine  fpx__^ 

a  fast  sort  package.  (Merge-ordering  uf  rwu  groups 'which 

are  in  themselves  ordered  can  be  e„asil*y  done  in  no  more  > 

than  M  operations.)  Thus  this  algorithm  should  be  Useful 

'  2  ' 

as  a  general  procedure  for  speeding  many  N  calculations l 
The  algorithm  was  applied  in  this  case  in  order  to 
avoid  certain  difficulties' mentioned'  in  the  next  section, 
not  out  of  considerations  of  speed.  Nevertheless,  it  is 
interesting  to  note  that  for  1-arge  N  a  more  efficient 
combining  rule  than  (2)  can  be  used.  With  the  natural  ^ 

identification  of  P(^|S)  as  the  elements^df  an  array, 

(2)  is,  formally  equivalent  to  a  convolution.  Well-known 
convolution  techniques  us,ing  the  Fast  Fourier  Transform 
could  be  uVpd  to  reduce  the  total  number  of  computations 
to  the  order  of  N(log2N)2.  Due  to  the  fact  that  the  FFT 
is  more  efficient  only'  for  rather  large  N,  coupled  with 
the  necessity  of  padding  to  avoid  wraparound?  t'his"  „ 
approach  was  not  thought  desirable  here.  For  extremely 
large  N,  it  might  hold  some  merit,  but  in  tha'-t/  cdse  v 
limiting  approximations  might  be  more  appropriate  any-  , 
way.  / 


J 
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COMPARISON  WITH  PREVIOUS  METHODS 


Feller  (1950)  gives  a  formal  solution  for  the  prob¬ 
ability  of  exactly  k  events 


l  N 


poo.  =  £  (-D3_k  (?)  e 


j=k 


K  3 


(3) 


where 


N 

*1  ■  S  ^ 

x  i=l 


E2  s  £  PiPj 

1^3 


(4)' 


etc. 


*  . 


Booker  (1965) ,  in  order  to  avoid  the  computation 


in  (4),  transformed  this  result  into  polynomials  in  Sj  •, 
where  ■  •“ 


N 


^  [Pi/u-pp] 


S. 


i 


with  coefficients  related  tp  the  Stirling  numbers.  His 
method  may  be- useful  for  calculating  only  a  few  of  the 
P(k),  for  small  k,  but  if  afl  the  P(k)  are  to  be  calcu¬ 
lated,  then  it'  will  be  quite  slow,  even , assuming  that 
all  £he  coefficients  are  either  available  or  may  be 
recursively  calculated  (not  an  easy  task  in  itself)'. 

In  addition,  the  method  will  not  work  very  well  for 
close  to  1,  due  to  the  division  in  (S) ,  nor  is  it 
obvious  that  his  alternating  series  are  numerically 
well -conditioned. 

This  author  (Wirth,  1970)  showed  that  all  of  the 

E.  could  be  computed  at  once  recursively  in  approxi- 
3 

mately  N  /2  operations,  using  the  scheme 


E1  Pl*  E2CP2E1'  E3  P3E2'  **"  EN  PNEN-1 

®l“®i+P2''  ^2_E2+P3E1'  E3_E3+P4E2*  ***'  ^N-l— EN-l"^PNEN-2 

ErEl+P3'  E2=E2+P4E1'  E3=E3+P5E2'  ***'  EN-2=EN-2+PNEN-3 


El=El+PN-r  E2"E2+PNE1 


ei  ei+?n 


(read'ljke  print,  left  to  right,  top  to  bottom).  With 

-  ■  *  / 

this  method',  all  .gf  the  P(k)  may  be  calculated.  $t  once, 

in  ~N  operations.  The  resulting  program  (NETPROB  I, 
Appendix)  is  remarkably  s imple ,  but  suffers  from  a  ”v 


f 

1 


major  limitation:  for  all  clo..e  to  1,  the  terms  in 
(3) '  can  become  very  large  and  truncation-'  effects,  can 
become  serious.  The  author  showed  that  the  largest  term 
will  be 

'v*1 

,  * 
r~~v 

N  i  Nl  N 

max  (“)  ( J)  «  - r  »  3*7  N 

•  [(N/3)!]3 


„  O 

For  N  =  20  this  is  ~2  x  10  and  the  worst  ca^e  accuracy 

f  '  '  - 

obtainable' with  CDC  1604  36-bit  floating-point  mantissas 

is  only  —1/2% . /Total  loss  of  significance  occurs  for 

N  =  27.  Use  of  double-precision  would  allow  the  range 

of  N  to  be  extended,  but  is  extremely  time-consuming. 

Nevertheless ,- for  small  N  the  prograiri  ls  eminently 

satisfactory.  (The  recursion  for  the  E.  should  also  be 

'  ^  { 

useful  for  reconstructing  a  polynomial  from  its  rootp. ) 

The  method  described*  in  the  first  part  of  the  pkper, 
however,  suffers  from  no  such  limitation  and  is,  in  addi¬ 
tion,  slightly  faster.  Thus  it  appears  to  be  distinctly 
■'  v~~\ 

preferable.  Tests  performed  by  comparing  NETPROB  III.  ■ 
against  a  double-precision  version  of  NETPROB  I  for  all 
N  <_  40  showed  excellent  agreement  and  no  loss  of  sig¬ 
nificance. 
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APPENDIX 

Computer  subroutines  are  written  in  FORTRAN-63,  a 
programming  language  of  the^CDC  1604B  computer. 


I 

S 


1 


SUrROL  T I Mfc  N  E  I  P  fi  0  B  (  N,P,  PN  >  III 

DIMENSION  P  ( N  )  ,  Pn<N>.  WI200> 

PN<K>  =  pRoh.  Up  AT  LEAST  *  OUT  Of  M  E  VtN  T  S  “HEN  INDIVIDUAL  EVENTS 
HAVE  FRO^.  P<1>.  N  5fl  ft.  wIRTH 

II  *  2 

DU  5  1  z  1,n 

h 1 1 1 >  =  p(i>  s  wni-i>  p  i p(d 

5  II  *  II  ♦  2  ' 

MNG  *  2  S  LOC  s  NT2  =  2  *  N 

IV  MNG2  -  MING  S  N2  =  MNu2  /  2 

MBS  *  NING  *  2  T  NNG  =  NT2  -  NlNG 

LUC2  s  LOC  $  LOC  =  N  f 2  -  LOC 

I F  <  NNG  )  100»100»2'i 
2t  OU  30  1 8  =  i , NNG »MNO 
L  =  LCC  +  IB 

.CALL  FHOrK(  N2*R<L)<L(l  +  NING2>*  NI NG2' * (L0C2+ IU  )  > 

3o  II  =  lb 

IB  =  U  ♦  NING  S  NR  =  NT2*1-IB 

1 F  (  Nh.GT.MNG2  )  4U.5U  _  _ _ 

4o  L  s  LCC  *  IB 

CALL  FKOBiU  I)f2»  L  <L  )  *  L  U.  +  MNG2  >  •  NH'2,  W  ( L0C2+  I R  >  > 

G(-  TO  in 

5U  Du  60  1  s  Ib,M2 
6v  K(LOCi+I  )  =  w«LCC*I  > 

Gl  TO  ip 

lOu  CALL  PRORM  n2/L(L0C*.1),W(L0C  +  N|NG2  +  1)»  N  »  W  <  L0C2*1 )  ) 

K  s  n  %  PN(K)  s  w l L0C2+K+1 ) 

DO  200  1  =  ?,N 
K  =  K  -  1 

200  F  N  <  K  )  =  ■PNIk  +  iJ  ♦  W‘L0C2+K+1) 

R6TUHN 

END 
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160 

170 

180 

190 

200 

210 

220 

230 

240 

25o 

260 

270 

28o 

290 

30  0 
31o  ’ 
320 
330 


SUpRoC I i Nfc  PPObK (  M2,P2,Pl*  N I » P T  ) 

DIMENSION  Pi<N2>,  P2lf'2>»  P  T  <  N  T ) 

FT  (K  +  i )  3.  PRnbi  Oh  tXACTLV  «  GUI  OF  NT  EVENTS  T  0  •  A  L  FROM  TWO 
GROUPS#  VITh  Pl#P2<l*l>  »  MROd-  OF  E*«CTlT  J  EVfcNTS  OUT  OF  Ni»N2* 

Nl  s  NT  -  Ng  $  IS  =  0 

FT  a  H  *  pj 
Of-  100  K  s  i. NT 
P  T ( *♦! J  1  0  . 

IF (  K-GT.Ni  )  20*l0 
10  il  *  n  \  Go  jo  40 

20  IL  s  Ni 

IF<  K.GT.N2  >  3C » 4y 
30  IS  !  K  -  N? 

40  DO  100  I  =  |S,  1L 

100  F 1  ( K  ♦  j  )  =  PT(K  +  i)  ♦  Fl<  I«l>*P2<K-In>  1 

RETURN 
END 
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l4o 

150 
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100 


SUBROUTINE  Ng  TPfiOU  <  N,P,  pN£T  ) 
DIMENSION  P  ( jg )  *  oliei<N),  g.(  2  0  ) 

PtogT IK )  a  PPOb •  OF  AT  LpaS*  K  OUT  UF 
EVENTS  MAVg  pnOB.  PlJ).  N  £  2t) 


N  events  when  individual 

M.  W I  rt  Tm 


l 


CALL  ERASE!  N»t  ) 

DO  10  K  a  1,N 

E  s  E  ♦  P(K)  $  jm  *  N  « 

DO  10  J  a  l.jM 
10  E(j  +  1)  =  £(j  +  l)  ♦  P<<’. 

PriET(N)  a  fc(N>  \  NM  8  N  " 

DO  30  K  a  l.NM 

KR  •  N  _  K  S  PNE I IKK) 

£  8  **  *  .c  ■  ” 0  • 

DO  20  J  a  KP,N 

D  8  B  +  i  •  t  c  a  C  ♦  i 

20  RNgTlKRJ  =  PNbUxR)  v  a*E  ( J) 

30  PNgT (KR>  a  pnBI(kH)  ♦  PNgT<KP*D 
RETURN 

end 


K 

1 
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